function compoundAP(varargin)


% This computes a compound action potential from summing individual action
% potentials based on the product of two Boltzmann equations to give a 
% realistic shape. The action potentials are given a random delay 
% (representing conduction velocity) based on a mean +/- SD described 
% by a gaussian distribution with a skew factor. This reflects the
% skewed distribution of axon diameters in the sciatic nerve.
% 
% 2000 of these action potentials are summed to create the compound AP. 
% Slowing is introduced by slowFac being >1.0. This slows the conduction
% velocity of every action potential in proportion to its delay.

% Requires the Statistics and Machine Learning Toolbox Add-On for
% pearsrnd() and pearspdf()


%default values

save=false;
mean=2;
stdev=0.5;
skew=0;

% allow user input on command line: skew, then SD, then save or not
% two files created for save, the raw data for the compound AP, and a
% parameter file with the rise time to 50%, fall time back to 50%, 10-90% 
% rise time, AP width at 50% and peak amplitude and peak time
% of the compound AP

switch nargin

    case 3
           if  isnumeric(varargin{1}) && varargin{1}>-1 && varargin{1}<1
            skew=varargin{1};
           else
            skew=0;
           end

           if  isnumeric(varargin{2}) && varargin{2}>=0
            stdev=varargin{2};
           else
            stdev=0.5;
           end 
    
           if varargin{3}==true || varargin{3}==false
                 save=varargin{3};
           else
            save=false;
           end
    case 2
            if  isnumeric(varargin{1}) && varargin{1}>-1 && varargin{1}<1
            skew=varargin{1};
           else
            skew=0;
           end

           if  isnumeric(varargin{2}) && varargin{2}>=0
            stdev=varargin{2};
           else
            stdev=0.5;
           end 
            
           save=false;

    case 1
           if  isnumeric(varargin{1}) && varargin{1}>-1 && varargin{1}<1
            skew=varargin{1};
           else
            skew=0;
           end
 
        stdev=0.5;
        save=false;


    case 0
        skew=0;
        stdev=0.5;
        save=false;
end

timeArray=[0:0.01:5];

rng(0,"twister");
num=size(timeArray,2);

% compute the delays based on mean +/-SD

delays=pearsrnd(mean,stdev,skew,3,2000,1);  % skewed normal distribution
delayDist=pearspdf(timeArray,mean,stdev,skew,3);  %kurtosis is kept at 3

myData=zeros(num,5);
myData(:,1)=timeArray;
slowFac=[1.0,1.05,1.1,1.2];

for i=1:4

    apCompoundWaveForm=zeros(1,num);
    for c=1:2000

         %compute the rising and falling phases using Boltmann equations
        
         risePhase=1./(1+exp((timeArray-delays(c)*slowFac(i))/-0.05));

         fallPhase=1-1./(1+exp((timeArray-delays(c)*slowFac(i)-0.3)/-0.2));

         %now multiply them to get a single AP
         apWaveForm=risePhase.*fallPhase;

         %add together to generate the compound AP
         apCompoundWaveForm=apCompoundWaveForm+apWaveForm;

    end

    myData(:,i+1)=apCompoundWaveForm;
end


subplot(5,1,1)
plot(timeArray,risePhase);
title("Rising phase");
ylabel("Amplitude");

subplot(5,1,2)
plot(timeArray,fallPhase);
title("Falling phase");
ylabel("Amplitude");

subplot(5,1,3)
plot(timeArray,apWaveForm);
title("Individual AP");
ylabel("Amplitude");

subplot(5,1,4)
plot(timeArray,delayDist);
titleStr="Distribution of delays: skew " + num2str(skew)+" SD " + num2str(stdev);
title(titleStr);
ylabel("Probability density");

subplot(5,1,5)
plot(myData(:,1),myData(:,2));
title("Compound AP");
hold on;

plot(myData(:,1),myData(:,3));

plot(myData(:,1),myData(:,4));

plot(myData(:,1),myData(:,5));
xlabel("Time (ms)")
ylabel("Amplitude");
hold off;

if save==true
    % saving the raw data
    filename="CompoundAP_skew_"+num2str(skew)+"_SD_"+num2str(stdev)+".csv";
    writematrix(myData,filename);
end

%now do some measurements on compound AP

riseT=[0,0,0,0];
fallT=[0,0,0,0];
peakTime=[0,0,0,0];
maxVal=[0,0,0,0];
amp10=0;
amp90=0;
time10=[0,0,0,0];
time90=[0,0,0,0];
time1090=[0,0,0,0];
    
for j=1:4
        
    [M,I]=max(myData(:,j+1));
    maxVal(j)=M;
    peakTime(j)=I;
    halfAmp=round(maxVal(j)/2);
    amp10=round(maxVal(j)*0.1);
    amp90=round(maxVal(j)*0.9);
    firstCross=false;
    amp10Cross=false;
    secondCross=false;
    amp90Cross=false;
    for i=1:num

       if myData(i,j+1)>amp10 && amp10Cross==false
            time10(j)=timeArray(i);
            amp10Cross=true;
       end

       if myData(i,j+1)>amp90 && amp90Cross==false
            time90(j)=timeArray(i);
            amp90Cross=true;
       end
   
        
       if myData(i,j+1)>halfAmp && firstCross==false
            riseT(j)=timeArray(i);
            firstCross=true;
       end
       if myData(i,j+1)<halfAmp && firstCross==true && secondCross==false
            fallT(j)=timeArray(i);
            secondCross=true;
       end 
    end
end
width=fallT-riseT;
time1090=time90-time10;

%output to command window

for i=1:4
    timings=[riseT(i),fallT(i),time1090(i),width(i),maxVal(i),timeArray(peakTime(i))];
    disp("slowing factor: "+num2str(slowFac(i)));
    disp("Rise time|Fall time|10-90 Rise|AP width|Peak amplitude|Peak time");
    disp(timings);
end

%optional output to file

if save==true
    % creating the parameter file with measurements
    filename="CompoundAP_skew_"+num2str(skew)+"_SD_"+num2str(stdev)+"_params.csv";
    for i=1:4
        timings=[riseT(i),fallT(i),time1090(i),width(i),maxVal(i),timeArray(peakTime(i))];
        lineStr="slowing factor: "+num2str(slowFac(i)); 
        if i==1
            writematrix(lineStr,filename);
        else
            writematrix(lineStr,filename,'WriteMode','append');
        end
        lineArray=["Rise time","Fall time","10-90 time","AP half width","Peak amplitude","Peak time"];
        writematrix(lineArray,filename,'WriteMode','append');
        writematrix(timings,filename,'WriteMode','append');
    end
end

end

